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TECHNICAL MEMORANDUM X-64658 


THE SELECTION OF APPROXIMATING 
FUNCTIONS FOR TABULATED NUMERICAL DATA 

INTRODUCTION 


In scientific work the scientist is often confronted with having to use 
tabulated numerical data. In many cases, it is desirable to represent the tab- 
ulated data by an approximate function containing arbitrary constant coeffi- 
cients, so as to obtain the "best” fit to the given data. This report is concerned 
with the development of a computer program that selects (from a list of can- 
didate functions) the approximating functions and associated coefficients which, 
when combined into a weighted sum, results in the "best" fit of the numerical 
data. 


The "best" fit of a set of numerical data is defined to be the one which 
minimizes the sum of the squares of the differences in the tabulated values and 
the corresponding values of the independent variable. Elementary matrix 
operations and vector methods are the techniques used to determine the recur- 
sion relations which yield the coefficients of the selected approximating func- 
tions. The advantages of this "curve fitting" program are as follows: (1) 
multivariable approximations can be performed, (2) flexibility with respect 
to the type of approximations used, (3) the program is designed to choose 
the "best" terms to be used in the approximation from an arbitrary list of 
possible terms so that little knowledge of the proper approximating form is 
required; and (4) recursion relations are used in determining the coefficients 
of the approximating functions, which reduces the computer execution time of 
the program. 

Included is a discussion of the problem, a method of solution, deriva- 
tion of the recursion relations of the coefficients, mathematical justification 
of the solution, a listing of the computer programs, and a section showing 
results from the use of the program. 


DISCUSSION OF THE PROBLEM 


From a set of n data points for m independent variables and one 

dependent variable, denoted as (y., x.,, x.„, . . . , x. ) where 

1 11 i2 im 

i = 1, 2, . . . , n, a table may be constructed similar to the following matrix: 




^12 

• • • 

^Im 

^2 

^21 

^22 

• • • 

^2m 

• 

• 

• 


• 

• 

• 

• 

• • • 

• 

• 

• 

• 


• 

y 

n 


X 

n2 

• • • 

X 

nm 


Such a table is assumed to define a function y - f (xj , X 2 , . . . , x^) 

over a region including the n data points. The desired approximation could 
be the sum of N specified functions where N is an arbitrary integer. In this 
case, the approximating function would have the form y^ = aj^f^ (x^, X 2 , . . x^ 

+ a 2 f 2 (xi , X 2 , x^) + ... + a^f^ (xi , X 2 x^) where, 

f . (xi , X 2 , . . . , x^) , for i = 1 , 2 , . . . , N are arbitrarily specified 

functions which must be defined over the region including the n data points. 

The coefficients, a. , of the specified functions in the approximations are to 

be determined as similar to those in "curve-fitting" problems. The first 
task accomplished by this report establishes the "best" L (where L is an 
integer such that 3 < L < N) approximating functions which are chosen from 
the total list of N possible functions. This allows much greater flexibility 
thfin the standard "curve-fitting" approaches because many different types of 
functions can be included in the list of N functions, and a selection of L func- 
tions, which produce the "best" fit to the data, is made automatically. In 
Reference 1 a similar problem is solved but the computational scheme devel- 
oped here is much more efficient than the one described in Reference 1. In 
the next section the conditions which determine the coefficients for a fixed 
number of approximating functions are derived. In the following section recur- 
sion relations are derived, and the method for selecting the "best" L terms 
is explained. The succeeding section contains the numerical results followed 
by the conclusion. 


METHOD OF SOLUTION 


Once the tables from the n data points of m independent variables 
and one dependent variable have been made and the choice of the N arbitrary 
functions, to be used in the approximation has been decided, a second table 
may be constructed similar to the following matrix: 
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*11 

'21 

• • • 

^N1 

^2 

'l2 

*22 

• • • 

^N2 

^3 

^13 

*23 

• • • 

• • • 

^N3 

^n 

*ln 

*2„ 

• • • 

• • « 

^Nn 


th 


til th 

where, the i'"** function (f ) evaluated at the j data point is denoted as f.. 

(1 < i < N and 1 < j < nj . Now, each of the columns in &e preceding table 
may be used to define a vector as shown below: 




"*ii’ 


’*21' 


^N1 

^2 


^12 


*22 


^N2 

^3 


*13 


^23 


^3 

« 

• 

^1 = 

• 

• 

^2 = 

• 

• • • " 

• 

^n 


_*!„_ 


>_ 


^Nn 


The set of vectors g. , where 1 < i < N , is assumed to be a set of linearly 

independent vectors. Thus, after the choice of the N arbitrary functions to be 
used in the approximation has been determined and the components of the set 
of vectors denoted as g^ have been calculated, an approximation to y may be 
written as 

y^ = aiii + aaiz + . . . + 


3 



or 


N 



a. g. 
1 1 


Now, define a quantity D(a) in the following manner: 
Let 

D(a) = (y-ya^ ‘ 


( 1 ) 


( 2 ) 


where, " denotes the dot (scalar) product of vectors and 


^2 

^3 




Note: D(a) as defined by equation (2) is an indication of the error between 
and the approximation y . By substituting equation ( 1) into the right member 

of equation ( 2) , D( a) may be written as 


D(a) = (y - J a. i.) . (y - 2 ^ sj » (3) 

i==l i=l 

and since the dot product of vectors is distributive over addition, equation ( 3) 
becomes 

N N N _ _ 

D (S) = y . y - 2 2 a (y . g ) + ^, ^ a a (g. . i ) . ( 4 ) 

i=l i=l j=l •’ 


4 


'< I 



Before continuing, matrix notation will be employed to simplify the 
expression involved. Equation (4) may be expressed in terms of matrices 
as will be shown. A colunan vector, for example a , may be considered as a 
matrix whose dimension is N x 1 . Similarly, a row vector may be considered 
as a matrix whose dimension is 1 x N . Let G be the N x N matrix whose 
elements are dot products of the vectors, 1 < i < N ; i. e. , 


] 


with i = 1 , 2 , 


N and j = 1, 2, 


N 


Furthermore, define the following quantities as shown below. Let 

y • ii 
y • g2 


_T 

and, with z denoting the transpose of the N x 1 matrix z ; i. e. 


j^y • ii y • §2 • • • y • 

and similarly 


-T 

Thus, by use of the above definitions, equation (4) may be written as 

— —X — “X “ —X “ 

D(a) = yy- 2 az+aGa . ( 5 ) 


5 


It is desired to minimize D(a) , which is an indication of the error between 

V and V with respect to the coefficients, a. (i = 1 , 2, . . . , N) , or with 
•' a 1 

respect to a. To show that this may be accomplished, D(a) will be expanded 
in a Taylor series about an arbitrary point, a* , where 


aj* 



The Taylor series will have the form 


D(a) = D(a*) + 


ao(a) 

0a 


(a - i*) 


a=a* 


+ -;(a - 


0a 


0D(a) 


0a 


(a - a* ) 


a^a* 


+ higher order terms 


( 6 ) 


The partial matrices in equation ( 6) can be readily computed from equation 
(5) . Therefore, by use of equation (5), 


and 


0D(a) 

0a 


-T -T 
-2z + 2a G 



( 7 ) 


( 8 ) 


6 


Thus, the higher order terms vanish. That is, all terms involving the 
partial derivatives of D(a) , where k = 3, 4, . , . , are zero because G is 
independent of a or equivalently, the a.(i - 1, 2, . . . , N) . Hence, the 

Taylor series of D(a) , equation (6) , may be written as 



(9) 


and by substitoting equations (7) and ( 8) into equation ( 9) , the Taylor 
series of D(a) becomes 


D(a) = D(a*) + 



+ 


2a* 



(a - a*) 


+ 


1 

2 



2G(a 


a* 


( 10 ) 


hi the following discussion it will be shown that equation (10) is 
equivalent to equation (5). That is, the truncated Taylor series of D( a) 
about the arbitrary point a* is a precise representation of D(a) . Thus, 
the minimization of equation (10) is equivalent to the minimization of D(a) . 
To show this, consider the right member of equation (10) which is shown 
below: 


D(a}=) + 


(- 


-T 

2z 


T 


■) 


+ 2a«< G/ (a - a*) 


+ 




( 11 ) 


When all the terms are expanded, expression (11) can be written as: 


7 


( 12 ) 


D(a*) 


-2z'^a + 2a*^Ga + 2z"^a=- -2a*^Ga* + a^Ga 
- a* ^Ga - a^Ga* + a* ^Ga* 


rn 

Since a* Ga may be considered as a 1x1 matrix (a scalar) it is equal 
to its transpose. Therefore 



but = a« and since G is a symmetric matrix 

- T - -T - 

a* Ga = a Ga* 


Thus, by collecting like terms, expression (12) becomes 


D(a*) - 2z^a + 2z^a* - a*^Ga* + a^Ga 


(13) 


Now, an expression for D(a=;=) may be obtained from equation (5) 
expression (13) may be written as shown below: 


Thus, 


_ _ _ 'p_ _ T> _ -T_ -T- - T - -T - 

y . y - 2a* z + a* Ga* - 2z a + 2z a* - a* Ga?^ + a Ga 


(14) 


T- 

As before, because a* z is a 1x1 matrix, it is equal to its transpose. 
Thus 




8 



since 



and by collecting like terms, expression (14) becomes 
-T- -T - 

y . y - 2z a + a Ga , 


which is D(al as shown in equation (5). Thus, equation (10) is equiva- 
lent to equation ( 5) . 

As mentioned previously, it is desirable to minimize D(a) [equation 
(5)] with respect to a; but since it has been shown that equation (10) is 
equivalent to equation (5), equation (10) will be minimized with respect 
to a . Recall that 


D(a) = D(S*) + 



(a - £*) 



(a - S*) 


where 


9D(a) 

3a 


-T -T 
-2z + 2a G 


and 



2G 


9 


To minimize D(a) as given by the above expression, consider the two cases 
where 


3D(a) 

. ^ . a=a* 

and 

8D(a) 

9a — 

L J a=a* 

In the first case where 



it can be seen that the vector a* will produce the minimum value of D(a) 
equal to D(a*) if it can be shown that the quadratic form 



is positive for any choice of a other than a = a* . To see tl^t this quadratic 
form is, in fact, positive for any choice of a such that a =^a* note that 



10 



as shown in equation ( 8) . Thus, 



-.vT 



’aD(a) 

\ 

1 da 

da 



(a 


a=a« 


a*) 


-T - 
X Gx 


(15) 


where x = a - a* ^ D, Now, the quadratic form x Gx, can be written in 
its matrix form as follows; 


— T - 

X Gx = ^Xi, X2, . 


•"n] I 




“ • 

Si • Si Si • S2 • • • Si • Sj^ 


Xl 

S2 * Si S2 • S2 • • • S2 • Sjj 


X 2 

• • • 
• • • • • • 


• 

• 

• • • 

Sj^ * Si Sjj • S2 • • • Sj^ • Sj^ 


• 

m m 


(16) 


where x^ = a^ - for i = 1 , 2 , . , . , N . By performing the indicated 
operations in equation (16) , the right hand side can be expressed as 


N • N 

Z X Ag. . ii) Z X • S2> • 
3=1 ^ ‘ 3=1 ^ ^ 



11 


which is equal to 


N N N _ _ 

2 X (i . ii) Xi + 2 X <i . iz) X2 + . . . + 2 X (i . i ) X . 

_j=i ^ ^ j=i ^ ^ j=i ^ ^ 

(18) 

Since the x.(j=l, 2,. . . ,N) are scalars, expression (18) may be written as 


N N N 

2 (XlSl • X.g.) + 2 (Xzg2 • X.g.) + . . . + 2 

j=l J J j=l ^ ^ 3=1 ^ ^ j j 


(19) 


and because the dot product of vectors is distributive over addition, expression 
(19) becomes 


( Z 

\j=‘ ' 7 


Vn> 



,( 20 ) 


which is the square of the magnitude of the vector. 


N 

Z 


i=i 


x.g. 

] ] 



that is 


N 


j=l 


x.g. 
J ] 


(21) 


12 



Thus, 


— T — 
X Gx 


N 


Z Yi 

j=l ^ ^ 


> 0 


f 


since for some i such that 1 < i ^ N and x. ^0 and the set of vectors 

1 

g. (i = 1 , 2 N) is a linearly independent set. Therefore 



a*) 



a=a* 


(a 


a*) 


is always positive for arbitrary a^^a* . So, when 


8D(a) 

da 


= 0 


a=a* 


the desired minimum of 


D(a) = D(a*) + 


9D(a) 

9a 


(a -a*) 


” ~5jC 

a=a 




_9_ 

9a 


9D(a) 


9a 


(a -a*) . 


- -9ii 

2e=a 


is obtained at a’‘‘ . 


13 


Now, to investigate the second case. 


9D(a) 

9a 


# 0 


a=a* 


( 22 ) 


it will be shown that, should 


9D(a) 

9a 


# 0 


a=a* 


9 


then D(a) cannot obtain a minimum at a* . This will be shown to be the case 
by proving the following. If 


9D(a) 

9a 


0 


a=a* 


9 


A 

then there exists a ^ a* such that 


9D(a) 


9a 


a=a* 


A 

(a 


a*) < 0 


and 



’9D(a) 

A 

(a - a*) 

> 

1 A - T 

j (a - a*) 

u 

9D(a) 

1 - 

> (a - a*) 


9a 

a=a={« 



( 9a 

9a 

ia.a 


14 


to |m- 


Once this result is proven, then by equation (9) it is easily seen that 


A 

D(a) = D(a*) + 



(a - a*) 


1 - - T 

2 (a -a") 



aD(a) 

9a 


A 

(a - a* ) 


a?=a 


A 

which_implies that D(a) < D(a*) ;_therefore, D(a) does not attain a minimiun 
at a=a* . To show this, suppose ai^a* and let g be the absolute value of the 
largest element in the N x N matrix. 


( 

_ ■ 

) 

) 9 

9D(a) 


( - 



i 9a 

9a 

) 


Let N (a*) be an €- neighborhood (c > 0) of a« and aeN (a*). There- 
€ _ _ _ _ e 

fore, I a - a* I < e but, a - a* ^ 0. Then 



( 


\ 

N^ge^ > 

1 -^,T ; 9 

_(a - a=^) /_ 

9D(a) 

> (a - a*) 


/9a 

9a 

1 - - 


\ 


/ a=a* 


(23) 


The equality holds only when g = 0 but, since it has been shown that the term 
involving 



is positive, the inequality is strictly ’’greater than. ” 


15 


To see that inequality (23) is true, note that the ri^t side is a 
scalar and g is maximiun of the absolute values of fee elements in 



dD(si) 

1 

\ ^ 

9a 

j 


a=a* 


Thus, by performing the indicated operations on the right side, a summation of 
terms results. The left side of inequality ( 23) may be written as the 
summation of terms as shown below; 


-N^ge^ = -(g€^ + ge^ + . . . + ge^) . (N^ terms) 
2 2 


Thus, term-wise, the left side of inequality (23) dominates the right side. 
Now, since 


9D(a) 

3a 

J 


^ 0 


a=a* 


9 


there exists some component, say the 


.th 
1 , 


such that 


9D(a) 

9a. 


^ 0 


a = a* 


If 


9D(a) 

9a. 


< 0 

a=a* 


) 


16 


then it is possible to make 


3D(a) 

da 


(a 


a=a* 


a*) < 0 


by letting a. = a * + |- and a. 

1 1 2 j 

Now, should 


a.* where i 
] 


3D (a) 
3a. 


> 0 


a=a* 


then it is possible to make 


3D(a) 

3a 


(a - a*) < 0 

a=a* 


by letting a. = a.* 
11 

Thus, if 



a. = a.* 
J i 


where i 


3D(a) 

3a. 

1 


< 0 


a=a* 


and j = 1, 2, . . . , N . 


i and j = 1, 2 , . . . , N . 


17 




Pl> 




Therefore 


9D(a) 

A 

(a - a*) 


9D(a) 

. ^ . 

a=a=5' 


. . 




9D(a) 




9a 




X 


Thus, for any e, such that 


9D(a) 

9a. 

1 J 


€ < 


3 r— 3 , 


2N^g 


then 


9D(a) 


9a. 

1 


a=a* 


1 >N2g€2 


But, from inequality (23) , 


N^ge^ 


> 

i d - 1 — 

*9D(a) ' 

i (a - a*) 


^ (aa 

9a 

* a==a* 




20 


to I to|i^ 


and since 


_ _ 




r 


9D(a) 

A 



9D(a) 



(a - a*) 

— 




9a 

_ 



9a 



a=a* 



1 

a=a* 


as shown in equation (24) „ it follows that 


8D(a) 


9a 


— — .X. 

a=a'‘' 


A 

(5 - a*) 


> 



( 



1 ^ 

9D(a) 





I 9a 

. 9a . 



a=a ' 


A 

(a - a* ) 


(25) 


Since, 


A 

D(£) = D(a*) + 


9D(a) 


da 


a=a 


(^- a*) 



u 

9D(a) 

1 



j 


9a 

. 9a . 

I 


(a - £*) 


a=a ' 


A 

the result just derived shows that D(a) - D(a?!') < 0 . Thus, a* is not the 
point at which D(a) is minimized. Therefore, if a-' is the point at which 
D(a) is minimized, it is necessary and sufficient that 


9D(a) 

9a 


= 0 


a=a* 


21 



The preceding derivation and others to follow use many results found in 
many books on vector and matrix analyses. References 2, 3, and 4 are 
suggested for those unfamilar with these results. 

When the conditions previously derived to minimize D(ai are satisified, 
they may be used to compute the coefficients as shown below. From equation 
(71 


8D(a) 

8a 



-T 

+ 2a G 


Thus, it can be seen that 
the following relation: 


8D(a) 

8a 


will be the zero vector when a satisfies 


-T 
a G 


-T 

z 


which is equivalent to 


(26) 


N 

Z 


j=i 


aj(g; 




y • Si 


(27) 


Thus, by solving equation (26) or (27) for a. where i 1,2 

a unique set of values for a may be determined which will make 
These values can be used as a* so that 


, . • . , N, 

8D(a) 

8a 


0. 


8D(a) 


8a 


a=a* 


will be the zero vector. Then D(a*) will be the minimum value of D(a) , 
since it has been shown that 


22 



1 - - T 

I (a - a«) 


(a - a*) 


_a_ 

da 



a=a* 


is always positive. From equation (4), it is known that 


N N N 

D(a) = y . y - 2 J a (y • g.) + 2 2 ^ a (i . i.) 

i=l i=l 3=1 ^ J ^ J 


By factoring, it is possible to obtain 


D(a) 


N 

y • y + 2 

i=l 



a.(g. 


i.) 


2(y . gj) 


. (28) 


Since the necessary conditions to minimize D(a) [equation (27)] are that 


N 

2 • £) = y • gj . ( 29 ) 

j=i i •> 


they may be substituted into equation (28) to obtain 


N 

D(a*) = y . y - 2 a *(y . g ) , (30) 

i=i ^ 


where the elements a* are the values of a which satisfy equation (29) . 
Equation (30) may be used to indicate the error involved in the approximation. 


23 



At this point it should be noted that, for a fixed value of N , equation 
(26) could be used to compute the coefficients of an N-term approximation 

by use of the inverse of the matrix denoted as G ; i.e. , a = z [G] where 
[G] denotes the inverse of the matrfac G . Note that this is for a fixed value 
of N which, in turn, places a limitation on the type of approximation, since 
N is the number of terms and corresponding to this is the dimension of G . 

So, for an N-term approximation, the above equation could be used to compute 

the a (i = 1 , 2 , . . . , N) ; however, any other approximation (involving Nq 
i 

terms where Ng < N) would require computing the inverse of the G Matrix 
each time N is changed. This is true because each change in the number of 
terms of the approximation causes a corresponding change in the dimension 
of the G-matrix. 

As the number (N) of terms in the approximating function increases, the 
computer execution time of an algorithm using equation ( 27) to determine 
values of a* for an N + 1 term approximating function increases. For 
this reason, recursion relations have been developed which allow values of 
to be computed for an N + 1 term approximation using a large part of 
the work done in obtaining the coefficients ^ for the N-term approximation. 
Then the error in the approximation as given by equation (30) and the 
preceding recursion relations which will be derived in detail in the section 
to follow are combined into an algorithm for determining the 'Taest" L terms 
from a list of N possible terms to be used in the approximation (where 
3 < L < N) . During each step of this algorithm, an indication of the error 
in the approximation may easily be calculated by use of equation (30) . The 
first time through the computation, each of the gj.(i = 1 , 2 , . . . , N) vectors 

are xised as ii . The particular g. that results in the minimum error is then 

saved or stored as ii . The next time through the computation, each of the 
remaining g. vectors are used as & . As done previously, the g. resulting 

in the minimum error is then saved as g 2 . The remaining N - 2 vectors 
are then used as gg and the procedure is repeated. The process is continued 
iin t.ii the desired number L of the N (where 3 < L N) possible terms 
or functions have been chosen. 

RECURSION RELATIONS FOR THE COEFFICIENTS 

To derive general recursion relations which allow the computation of 

an n-term approximation to a given function (where n = 1 , 2 N) , the 

development for the cases n = 1 , 2 , and 3 will be shown in detail. Then 
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usmg the results for n = 1 , 2 , and 3 and also the results for n = 4 and 5 
which are not included, due to length and complexity, the general recursion 
relations which are listed can be obtained. 

To start this development, consider the case where n = 1 . Then 
equation (27) would have the form 


(gi • giXaii) = y • gi 
Therefore, the coefficient would be given by 

ail = , 

gl • gl Bii 


(31) 


where = gi . gj . 


n = 2 

For the two -term approximation, equation (27) has the following 

form: 




- 


■ - 

gl • gl gl • g2 


a2i 


y • ii 

g2 • gl g2 • g2 


a22 


y • g 2 


As before, it is desired to determine the coefficients in terms of g^ . gj , 
gl • g 2 > §2 • gl » §2 • §2 > y • ii > y • §2 • This is to be accomplished by use of 
the augmented matrix from the above system of equations. The augmented 
matrix is of the form; 


• 


- 

gl • gl 

ii • g2 

y • ii 

_g2 • ii 

g2 • g2 

y • g 2 
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where elementary row operations will be performed on the matrix to obtain 
an equivalent matrix of the system which will eventually yield the desired 
relations for the coefficients. The operations will be indicated to_the right 
of the matrix on which they are performed. Note that = Si • Si was 
computed in the previous section. From expression (32) , multiply row one 

of the augmented matrix by to obtain an equivalent matrix as shown below: 

Bii 


Si • ii 

gl . g2 

y . gl 


1 

Bi2 

an 

_ S 2 • Si 

g2 • §2 

y • g 2 _ 

Bn 

_ g2 • ii 

g2 • S 2 

y • g 2 . 


where denotes "is row equivalent to" and 


Bi2 


ii • g2 

Bii 


Note that a^ = ^ was computed in the previous section. Define 
Bn 

Bjj = §2 . ii and perform the elementary row operations as indicated below: 


■ 1 

Bi2 

an 


■ 1 

Bi2 

an 

. B 21 

g2 . g2 

y . g2 . 

^ /*D 1 

0 

B 22 

y . g 2 - a^ B 21 


-■r2-(B2i) rj 


where 

B22 = S2 ‘ S2 ~ B21 Bj2 
Finally, 


‘ 1 

Bi 2 

an 

/s/ 

■ 1 

Bi 2 

an 


B22 

y • g2 ■ ail B2 i 


f\ 

i 

yig 2 - an B21 

0 


0 

B22 


VB 22 / 


f 


which is the desired form of the augmented matrix. Hence, coefficients for 
the two -term approximation are given by 
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and 


a-22 “ 


y ♦ g2 - an Bu 
B 22 


> 


^21 - ^11 “ ®12 3-22 


n= 3 


Now, for a three-term approximation, the system of equations 
expressed in matrix notation is shown below : 


ii • Si ii • g2 

ii • S 3 ' 


3-31 


1 

1 

1 

§2 'Si S 2 • S 2 

S 2 • gs 


3-32 

= 

y • S 2 

. g3 • gl is • S2 

g3 • g3_ 


^33 


y ' S 3 


The augmented matrix may be written as 


gl • gl 

Si ' S2 

Si ' S 3 

y • ii 

S2 ' Si 

S2 ' S2 

S2 • S3 

y • g 2 

. S3 ' Si 

S3 ' S2 

S3 • S3 

y • g3 _ 


(33) 


Elementary row operations are to be performed on the augmented matrix, 
shown on line (33) , until it is in a form which will allow the coefficients asj, 
a 32 , and ass determined easily. As before, Bji = g^ • g^. Therefore 


gl • gl ii • g 2 ii • S 3 

y • il 


1 Bj2 Bis 

ail 

S 2 • ii g 2 • g 2 g 2 • is 

y . S 2 


B21 g2 . i2 g2 • S3 

y • S2 



- 


_ _ 

. g 3 • il S3 ' S2 S3 • S 3 

y • S 3 j 


[_ B31 gs . g2 gs . gs 

y • g 3 _ 
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where 


Bi2 


• S2 

Bii 


Bi 3 


Si > S3 

Bii 


an 


y • gi 

Bn 


B21 = g 2 • Si 
B31 = is • Si 

Note that B^2, B21, and a^^ have been computed previously; thus, only 6^3 
and 63^ are required to be computed here. Therefore 


'1 

Bi2 

Bi3 

an 

B21 

g 2 • g 2 

g2 • is 

y • is 

_B 31 

g 3 • g 2 

g 3 • gs 

y . is_ 


’ 1 

Bi2 

Bis 

an 

0 

B22 

g2 • gs “ B21 Bj3 

y . is - B21 ail 

_ B31 


_ 


is • g2 

gs • gs 

y • g3 
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with B22 - g2 • g2 - ^12^21 which has been previously computed. 


■ 1 

Bi 2 

Bi 3 

ail 

0 

B22 

S2 • S3 “ B21 Bi3 

y . g2 - B21 ail 

- B31 

S 3 • S 2 

S3 • Ss 

y • S3 


rs-Bsirj 


1 

Bi 2 


Bi 3 



^11 

0 

B22 

S 2 • 

S3 “ B21 Bi 3 

y . 

S 2 

- B21 ail 

0 

S3 • S2 “ B31 Bi2 

S 3 • 

S3 ■ B31 Bi 3 

y • 

S 3 

■ B31 ail 


Continuing, 

1 ®12 

0 B22 

0 §3 • §2 “ B31 Bi2 


Bi 3 

g2 • S 3 “ ^21 Bi3 
S3 • S3 “ B3I Bi3 


a^i 

y • S2 - B21 ail 
y • S3 " B31 ail 




^2 


■ 1 

Bi2 

Bi 3 

ail 

0 

1 

B23 

^22 

. 0 

B32 

S3 • S3 “ B31 Bi3 

y • S 3 - B31 ail . 


y • g? - B91 an 

with a22 = — , which has been previously computed, and defining 
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and 


„ _ g? • ga - B21 Bi 3 

^23 - 


B32 - 

is • g2 

- B31 Bi2 

Therefore, 

1 

Bi 2 

Bi 3 

ail 

0 

1 

B23 

as 2 

0 

B32 

is • is " ®31 Bi 3 

y . is - B31 ail 


r3“B32r2 


1 

B12. 

Bis 

an 

0 

1 

B23 

a22 

0 

0 

®33 

y • is - B31 an - B32 a22 


with 


B33 = is • §3 - B3I Bi3 - B32 B23 


Thus, 
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■ 1 

Bi2 

*Bi3 


0 

1 

®23 

^22 

. 0 

0 

B 33 

y • 6 ~ B 3 I ^11 “ ®32 ^22 - 



which is the desired form of the augmented matrix and where 



Hence, the remaining coefficients may readily be computed as shown below: 


^32 - 3^22 - ®23 ^33 


and 


^31 “ ^11 ■ Bi2 ^32 ~ Bi3 ^33 

The relations for four and five-term approximations were developed in 
a similar manner. Due to the complexity and length of the expressions 
involved, the development will not be given here. However, by reviewing the 
previous cases a recursion relation may be written as follows; for n = 1, 
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Bji - gi . gi 


ail 


^ y • g-l 

Bii 


For n 2 

B21 = S2 ‘ Si 


Bi2 


B 21 

Bii 


B22 - g 2 • S2 ~ B ^2 B21 


az2 


y ■ g, - ag B 2 ^ 
B 22 


a2i - a^i - a22 Bj2 

For i = 3,4, 5, . . .,N, where N is a positive integer, the relations are 
given by 


B. = g . g 
il ^i 


B 


B 


li B. 


U 

j 

11 
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B.. = g . g 
1 j 


B.. 



j-1 

2 


k=l 



B 

ik 


j = 2, 3, 


. , i-1 


B.. 

11 






a.. 

11 


i-1 


y • - Z a. . B 


k=l 


kk ik 


B 

ii 


a . . = a 
1 . i-J i-j , i-j 


Z B. . . . ^ 

1, i-j+k i-j, i-j+k 


j~ l«2,a«.,i-l 


As previously mentioned, and as indicated above, quantities used in determining 
the coefficients of a (k - 1) -term approximation are also used in determining 
the coefficients of a k-term approximation, (2 ^ k ^ N) ; hence, they need not 
be recomputed. 


To illustrate the use of the general recursion relations, let i = 3 and 
substitute it into the general relations to see if they do agree with the quanti- 
ties as calculated using elementary row operations on the three term matrix. 
For i = 3 , B^^ ~ given by the recursion relations. Therefore, 

Bsi" gs • gi > which agrees with Bgj as determined by elementary row opera- 
tions. Also, 
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hi 

Bll 


» 


by the recursion relations; therefore 


Bi 3 = but, Bgi = g3 • gi 


Therefore 


B 


13 


ga -..JS! 

Bll 


» 


which agrees with the results obtained using elementary row operations, 
continuing through the recursion relations, where j = 2,3,...,i-l 


B.. 

1] 



j-1 

S 


k?=l 



ik 


Therefore 

B32 = is • ia - ®12 ®31 


agrees with the derived value 


B.. 



9 


where j = 2,3,...,i-l. 


Then 
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Therefore 


B 


23 


^ 3 ? 

®22 


But 


®32 ~ S3 • §2 - Bi2 B 3I » 


therefore 


g = Sa^ « g 2 - B |9 B 3 ^ 

2d O 

^22 


where 


B31 = S3 . gi 


Bi 2 


gl • S 9 
Bu 


and 


Si . g2 = g2 • gl 
Therefore 


S3 . g2 


B 


fef) 


(S 3 • Si) g3 • S 2 - (g 2 • Si) 


23 


B 


i i A 
\ B^^ I 


22 


B 


'22 



and since, B21 - §2 * Si 



But 


i;i • il = il • Si and g:3 . g-2 = g2 . i:? 


Therefore 




S2 • Si 



B22 


where 


B 




gi ♦ Si 

Bii 


Therefore 


B23 " 


_ g? • Si - B,i B 


12 - 


B 


22 


which agrees with the derived value. 


B.. = 
11 



j-1 

2 \i=ik 
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Therefore 


®33 S3 . gs - Bj 3 B3J - B23 B32 


which agrees with the derived value. 


a 

u 


y 


^i 


j-1 

^kk ®ik 


B.. 

11 


Therefore 


a - y • S^ - ^ti - a.yy B^, 

^33 - R 

^33 


which agrees with derived value. 

The remaining coefficients are obtained from the following relation; 


a. . . = a. . . . - ) a . B 
1,1-J i-J.i-3 i,i-j+-k i-j,i-jt-k 

where j = 1, 2, ... i - 1. 

Therefore, a32 == a22 - B23 a33 and a3j == aj^ - Bj2 a32 - B13 a33 which 
agrees with the derived value. 

Thus, the recursion relations are valid for a three-term approximation, 
as shown above. If interested, one may derive the quantities for a four and 
five-term approximation in a manner similar to the one, two, and three-term 
approximations, which have been shown. As a result the reader can verify 
that the recursion relations hold for n = 4 and n = 5 . And an induction type 
of proof could be constructed to show that the recursion relations are, in 
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fact, true for any N > 3 where N is an integer. The Appendix is a listing of a 
printout of the programming of the recursion relations and the logic for select- . 
ing the "best" L terms from a possible list of N terms where (3 ^ L ^ N) . 

NUMERICAL RESULTS 


In order to show how the computer program works, examples of two 
types of approximations to a function of one independent variable will be in- 
cluded in this section. The numerical data to be used for these examples are 
a tabular function where tj(ETA) is the dependent variable and Mach is the 
independent variable. These numerical data are a particular set of aerodynamic 
data which are often used in trajectory simulation programs. In Tables 1 and 
2 the actual numbers indicating how 17 varies with Mach are given. In Figures 


1 through 20 , plots of this tabular data are given. For clarity, the following 
definitions are used to assist in reading Tables 1 and 2 : 

Y 

nominal value of the dependent variable 

XI 

value of the independent variable 

CY 

computed value as calculated by the approximation 

DIFF 

Y-CY; i.e. , nominal value minus the computed value 

PDIF 

percent difference, DIFF/Y 

DENOM 

value of the denominator; for polynomial approximations, 
DENOM = 1 

DDD 

the sum of the squares of the differences 

RMS 

root mean square; i.e. , DDD/M where M is the 

number of data points 

ERRl 

DDD/YDY where YDY is the sum of the squares of Y 

ERROR 

n/ ERRl/M 
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The two types of approximations to be considered are polynomial approx- 
imations (case 1) and rational function approximations (case 2). In each, 
a table printout and plotting of the approximations were initiated after the sixth- 
term approximation was made (i. e., NP - 6). For case 1, the sample print- 
out of all the polynomial approximations will be referred to as Table 1 and 
the sample printout for case 2 , the rational function approximations, will be 
referred to as Table 2 . The process was terminated after the fifteenth- term 
had been chosen by the program (NT = 15). 

Consider case 1 , that is, the polynomial approximations. The functions 
chosen to be used as terms in the approximations are listed below and may be 
seen in the program listed as the Appendix. 

G(l) = 1 

G(2) = X, where, x is the independent variable 

G(3) = -x^ 

G(4) = x^ 

« 

« 

G(19) = x^8 

G(20l = xl9 

So, from the above functions, the program will determine polynomial 
approximations ranging from one-term through NT-terms, as may be seen in 
Table 1 , where (NT = 15). The approximations will be of the form 


NT 

y = F(x) = 2 CA(J) . AFN(J) , 
J=1 


where AFN denotes approximating function number and CA is the computed 
coefficients of the corresponding approximating function. As may be seen 
from Table 1 , the following functions were chosen in the order shown below 
for a six-term polynomial approximation. 
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TABLE 1. POLYNOMIAL APPROXIMATIONS FOR 



OD^P • ,H5393HI5027*«6907-0I RMS ■ 1 2622559 1 282&>«<(92-.0 1 HM ■ , 1 5709 1 8875994339 + 03 


TABLE 1. (Continued) 
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TABLE 1. (Continued) 
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table 1 . (Continued) 
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TABLE 1. (Continued) 
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TABLE 1. (Continued) 
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TABLE 1. (Continued) 
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.1800*01 XI • .2800*02 CY • .17708882*01 BlfF- .29181878*01 PDIF« .18195378.01 DENOM ■ .10000000*01 
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DIFF* .11035977*00 POIF* .20838998*00 DENOM ■ .10000000*01 



TABLE 1. (Continued) 
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TABLE 1. (Continued) 
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TABLE 1. (Continued) 
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TABLE 1. (Continued) 
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TABLE 1. (Continued) 



69 




TABLE 1. (Continued) 




o o o ' o o 

o o o o o 

♦ ♦ ♦ ♦ ♦ 

O O o O k/t 

r> ^ 
» S' iTi kT 
if) if> • m ' io 4T 


Cl 

o 

<f 

o 

o 

cc 

un 


0 


o 

o 

♦ 

o 

o 

(N 

«G 


O 

o 

♦ 

o 

o 

o 

r*- 


o 

o 

♦ 

o 

o 

fy 

a> 


o 

♦ 

o 

o 

o 


o 

♦ 

o 

tv 


o 

♦ 

o 

w 

m 



o o o o 

♦ ♦ ♦ ♦ 

o o o o 

ty >0 O' «s 

if] Ip 


o 


o 

O' 

o 


o 

o 

< 4 ) 


o 

♦ 

o 

O' 

•€ 


o 

■* 

IT 

C 

rv 


X 


>■ >- >- V ► 


71 



6860192*01 POIFb . 9689765M*00 OENOM « •10000000*01 



72 



TABLE 1. (Concluded) 



o 9" in o — 

i/* r«fc f*v T' o o 

■o in a- fv o a 

r •m o «0 ' o O 

^ 0| IV IS. O O 

rv o ift o o 

o in *o M o o 

o rv >c j) o o 

a — n a o o 

o -• O O O O 

o a- 0 3 o o 

o u> III 7 a o 





ui 


o 

o 


o 


K Ul • (*> 
Z t I- FV 
o 

l«l »- X rn 
J >4 

c c 3 
■* o 

tK £ X M 
O X 
O X X X 
X «.X ol 
a. z uj 


o rv 

O IV 

m <n- 
in lii 
a- ^ 
o 

o o 

O (*> 

o o 

O 3 

o o 
o <« 


O 'O 
m — 
•O 1*1 

a* 

-• o 
< 

a r*> 

o >« 

O (V 

o -o 

a 7 

O Jl 


a iM 
o -e 
3 N 

a « 

o rv 

O ’C 

o <n 
a rv 
o in 
a rv 

O fv 
O -> 


rv o 
N O 
IV — 
fv N 
tv o 
O 

in n 
M — 
n in 

«*i 

z -• 
rg — 


o in 

a a 

o -I 

a a 
o o 
a a 

a o 
o o 
o o 
a a 
o o 
O o 


o a- 
•c o 
a 1*1 
'« o 
a a 
'O o 

a a 
•« a 
a o 
•« o 
a a 
'0 a 


— N tv tv O o 

■O y Jl .*v NO 

N O N NO 

“•1*1 O CM NO 

ON J> fM NO 

O m -O N NO 

a a N N NO 

o in T’ N NO 

On -on no 

ON [•) N NO 

ON ON NO 

O -• M -« NO 


N o 
a* <0 
-• o 

N O- 
<n o 
o o 


/*) T 
N M 
in 'O 
— T 
N O 
N « 


5 T ^ 
N O 

n o 
o 
u o 


o »• in < ^ o 

o »• NO TO 

O “• — • •• o 

o o -• o o 

a a T N a o 

O O IJI T O O 


r*> r) 
'© n 

CM 

Jl 

O -• 

o o 
o o 
o o 
a o 
o 

o o 


N N 
N — 
.'*» in 

j> <n 

N 


a in 
m >o 
a **> 
a N 

CJ N 
« — 


■o « 
n o 
Jl a 

N O 

M a 

— a 
-o o 
in a 
T a 
M a 

T O 


>1 a Jl O Jl -• 

a N '0(^0 »• 

O Jl <*l « O ^ T- 

fH T ^ 1*1 O T Jl 
•n 0-0 0 -ON 

a -• 04 o oi 


a a (*> Jl a 

o o 1*1 a 
o a rv »*» o 

a a -o Jl a 

a o T — o 

a o Jl o o 


.*1 o 

TOO 
N O UJ 
N a *- 
N o X 
-• o z 


r 

a iM N ox 

t*» r*> N O ;*! 

■n r o o »- 

■n f*» T o o 

N T o o yi 

N in t*i o Q 

Jl 

o a V o o >- 

a Jl :A a o- a 

a .N a a. o - -e 

a N o o o z 

•a ,N a a o < 

a •— -N a o: 



73 


RU(,IDi UA«PUT account: 3It62U PKOJICTI AVARnblNHli 

b337 - MAX TiMt 


AFN(l) is G(l) = 1 
AFN(2) is G(2) = x 
AFN(3) is G(3) = 

AFN(4) is G(4) = x^ 

AFN(5) is G(5) = x'* 

AFN(6) is G(10) = x® 

Thus, a six-term polynomial approximation (Table 1 and Figure 1 ) would 

have the form 


y = 


F(x) 


6 


= .E 

j=i 


CA(J) • AFN(J) 


or 

y = CA(1) = CA(2) X + CA(3) x® + CA(4) x® + CA(5) x^ + CA(6) x® 

where the coefficients obtained from Table 1 are given below. 

CA(1) = 0.3633551733860099 .+ 00 
CA(2) = 0.4452805637032931 + 00 
CA(3) = 0.4826696357426569 - 01 
CA(4) = 0.22136964 0 54 55174 - 02 
CA(5) = 0.3696033107572816 - 04 
CA(6) = 0 . 8652009407108565 - 13 
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If desired, the reader may determine approximations through 15-terms from 
Table 1 by following a procedure similar to the one just described for a 6-term 
approximation. Indications of the quality of an approximation are the graphs 
(Figures 1 through 20), and the error-term, as given on the printout. Table 
1 , for each of the 15 different polynomial approximating functions. The graphs 
of all the polynomial approximating functions follow and are identified as 
Figures 1 through 10. Figure 1 is a plot of the actual variation of tj with 
respect to Mach, which is shown as a solid line, and the 6-term approximation, 
which is plotted with asterisks. Figures 2 through 10 are similar plots for 
the 7- through 15-term approximation. 

As may be seen from Figures 7, 8, and 9, the 12, 13, and 14 -term 
polynomial approximations appear to be better approximations than the 6-term 
polynomial approximation. Hence, from Table 1 , it may be seen that the 
error does decrease as the number of terms in the approximations increases to 
14, at which time the smallest value of the error is reached (O. 5 x 10 ^) . Thus, 
of all the polynomial approximations considered in this example, the 14-term 
approximation would be considered the most accurate. From Figure 10 , it 
can be seen that the graph of the 15-term approximation consists of many 
"spikes." Also, from Table 1 (15-term approximation; AFN = 15) , note that 
the value of the error is approximately 0. 1033 . This is an increase compared 
to the 14-term approximation error and is due to a numerical precision problem. 
To overcome this, there has been developed a version of this program in which 
calculations are performed to 40 places [ 5] . 

Now, consider case 2, the rational function approximations to the same 
data considered in case 1. For this case, there are twenty functions designated 
for use in the numerator of the approximation and nineteen designated for use 
in the denominator. All of these functions are listed together below. 

G(l) = 1 

G(2) = X 

G(3) = x2 

G(4) = x3 

G(20) = x‘9 
G(21) = -xy 
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Figure 1. The six-term polynomial approximation. 
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Figure 2. The seven -term polynomial approximation. 
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Figure 3. The eight-term polynomial approximation. 
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Figure 6 . The 1 1 -term polynomial approximation . 
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Figure 7. The 12-term polynomial approjcimation . 
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Figure 8. The 13-term polynomial approximation. 
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Figure 9. The 14-term polynomial approjdmation. 
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Figure 10. The 15-term polynomial approximation. 
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G(22) = -xV 
G(23) = -xV 

G(38) = -x^V 
G(39) = -x‘V 


An explanation of how a part of the last 19 of the above functions became 
the denominator of a rational function approximation is given in the discussion 
to follow. 

The computer program listed in the Appendix determines rational func- 
tion approximations ranging from one-term through "NT" terms (where 
NT = 151 from the preceding list of 39 possible terms. A sample printout 
from this program is included here as Table 2 . Note that the symbols used 
in Table 2 are the same as those in Table 1 . As may be seen from Table 2 
the following functions were chosen (in the order shown belowl for a 10-term 
rational function approximation; 

AFN(l) is G(l) = 1 
AFN(2) is G(21) = -xy 
AFN(3) is G(22) = -xV 
AFN(4) m G(4) = x=^ 

AFN(5) is G(23) = -xV 
AFN(6) is G(5) = x^ 

AFN(7) is G(24) = -xV 
AFN(8) is G(6) = x® 

AFN(9) is G(25) = -x^ 

AFN(IO) is G(2) = x 
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TABLE 2. RATIONAL FUNCTION APPROXIMATIONS FOR 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 
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TABLE 2. (Continued) 


O 



114 



TABLE 2. (Continued) 
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TABLE 2. (Continued 
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TABLE 2. (Continued) 
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TABLE 2. (Concluded) 
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Hence, the approximation may be expressed as 'a function of the independent 
variable in the following manner: 

10 

y = F(x) = Yj CA(J) AFN(J) 
k=l 


where the coefficients are listed in Table 2. Thus, 
y - CA(l) + CA(2) • (-xy) + CA(3) • (-x^) + CA(4) . x^ 

+ CA(5) . (-x3y) + CA(6) . x^ + CA( 7) . (-xV) + CA(8) . x» 

+ CA(9) . (-x^y) + CA(IO) . x 


Therefore, the rational function approximation to the curve may be written as 


F(x) = V = + CA(lO)x + CA(4)x^ + CA(6)x^ + CA(8)x^ 

1 + CA(2)x + CA(3 )x2 + CA(5)x 3 + CA(7)x^ + CA(9)x® 


Again, as with the polynomial approximations, the plots resulting from the 
functions given in Table 2 are included here as Figures 11 through 20 . 

Note that the computed values are plotted as asterisks and the nominal 
values as dots. As may be seen by inspection of Figure 15 , the nominal 
curve and the curve as determined by the approximation almost coincide. Also 
note from Table 2 that the error for this 10-term approximation, 0. 66 x 10“^ , 
and that the algebraic sign of the denominator (DENOM) in the rational func- 
tion is always positive so that there are no zeros in the region of interest. 

Thus, this 10-term approximating fimction is considered the most acceptable 
of the rational function approximations. 

As before, on some of the graphs of the rational function approximations 
there are ’’spikes." Some of these are due to numerical precision problems, 
as is the case with the polynomial approximations, but not all of them. For ' 
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Figure 11. The six-term rational function approximation. 
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Figure 12. The seven-term rational fimction approximation. 
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Figure 13. The eigjit-term rational function approximation. 
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eta versus MACH 



ETA VERSUS MACH 


1.8 



Figure 15. The 10-term rational function approximation. 
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Figure 16. The 11-term rational function approximation. 
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Figure 17. The 12 -term rational function approximation. 
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Figure 18. The 13 -term rational function approximation. 
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Figure 19. The 14 -term rational function approximation. 
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Figure 20. Tlie 15-term rational function approximation. 
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example, see Figure 16 at xl = 10.5, From Table 2 it may be seen that 
the algebraic sign of the denominator changes from positive at xl 10.0 to 
negative at xl = 10. 5. Hence, there exists a zero of the denominator for this 
particular rational function approximation between xl = 10. 0 and xl = 10. 5 . 

As another example, from Table 2 , an error of 0.408 x 10“'^ occurs for a 
13-term rational function approximation [i. e. , at AFN (13)] which is 
less than the error of the 10-term rational function approximation. But, 
from Table 2 , it may be seen that the algebraic sign of the denominator is 
positive at xl= 6. 5 and negative at xl = 7.0 , indicating that a zero of the 
denominator exists between 6.5 and 7.0 for this approximation. Thus 
neither of these approximating functions can be used with any safety, although 
for the 13-term function both the error-term and the plot seem to be acceptable. 
Many other tabular functions have been approximated with the program. The 
sample printout and Figures 1 through 20 , should serve to familiarize the 
reader with the operation of the program. 


CONCLUSION 


The program, as designed, allows the user greater flexibility as com- 
pared to the standard "least-squares" approximation techniques. By use of the 
recursion relations to compute the coefficients of an approximation, quantities 
once computed need not be recomputed when the approximation is changed. 

Thus, many approximations to a function are easily and rapidly available. Also 
by use of the so called "picker" or "best g. choice," it is not necessary for 

the user to have a good idea as to the type of approximation that is desirable to 
approximate a function. The program chooses the best 1,2,..., N-term 
approximations to a function. Additional flexibility is available since the 
function used as terms in the approximations may be arbitrarily chosen. As 
shown by previous examples, pol 3 momial and rational function approximations 
are available for the same function during one run of the program. This 
was possible by setting ND = 0 and ND = 19 respectively in the data pack. In 
addition, the program is designed to accommodate multiple cases. 

The main restrictions are that the functions used as terms, to be chosen 
for the approximations, be defined over a region including the data points and 
that computer capacity not be exceeded by making the possible number of terms 
to be selected too large. 
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APPENDIX 

COMPUTER PROGRAM FOR DETERMINING COEFFICIENTS 
OF APPROXIMATING FUNCTIONS FOR 
TABULATED NUMERICAL DATA 
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Definition of Symbols 


Symbol 

CA 

GDG 

XDG 

G 

X " 

Y 

NC 

IPOINT 

NN 

ND 

NT 

TT 

NP 


Definition 

Coefficients of the approximation 

Dot product of the g, vectors 

Dot product of y and the vectors 

Arbitrary functions used in the approximation 

Independent variables 

Dependent variable 

Case Number 

Number of data points 

Maximum number of functions to be used in the numerator of 
the approximation 

Maximum number of functions to be used in the denominator of 
the approximation 

Maximum number of terms to be selected for the approxima- 
tion 

Constant added to the dependent variable 
Initiation of plotting 
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(JOlOO !• c multivariarle function approxination rith linear combinations or 
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A 

2 

00103 

5* 


IMPLICIT REAL48(A-H»0»Z) 


A 

3 

0010«l 

64 


REAL LPLOTX*LPLOTY»LPLOT2 


A 

4 

00105 

7* 


DIMENSION LABELXU21. LABCLYUZli LPL0TX(200 )i LPL0TY(200), 

1 LPLOTZ 

A 

5 

00105 

8* 


1(2001 


A 

' 6 

00106 

9* 


DATA /LABELX/6H MACH»6H TT ■ »246H |6H NO • 

,246H 

A 

7 

0O1O6 

104 


1 .ANKBL • »366H / 


A 

6 

00110 

1 !• 


DATA /LA6ELY/1246H / 


A 

9 

00112 

124 


CALL IDENT (9) 


A 

10 

00113 

134 


NC«1 


A 

11 

001 l«t 

144 


10 RCAD(St390,EN0a375»ERR*3751 1 PO 1 NT « PO I NT » NT > NN • NO , TT • NP 


A 

12 

00125 

154 


read (5,400) LABELY 


A 

13 

00133 

164 


ENCODE (520 iLABELX (41 ITT 


A 

14 

00136 

174 


ENCODE (530,LABELX (71 INO 


A 

15 

OOIMI 

184 


N«NN4N0 


A 

16 

00142 

194 


PRINT 380 


A 

17 

00144 

204 


PRINT 500, NC 


A 

18 

00147 

214 


PRINT 510, IP0INT,NN,N0,NT,TT,NP 


A 

19 

00157 

224 


DO 30 1«I ,N 


A 

20 

00162 

234 


DO 20 Jal ,N 


A 

21 

00165 

244 


20 G06( 1 , J)a0«0 


A 

22 

00167 

254 


30 Y06( I )a0«0 


A 

23 

00171 

264 


YOY«0,0 . 


A 

24 

00172 

274 


DO 100 Km \ , (POINT 


A 

25 

00175 

284 


READ 410, Y(K) ,Xn(,l 1 


A 

26 

00201 

294 


Y(»()*Y(X)4TT 


A 

27 

00201 

304 

C 

THESE ARE THE ARBITRARY FUNCTIONS CHOSEN 


A 

28 

00202 

314 


6( 1 ,K )■! • 


A 

29 

00203 

324 


6(2,K)>X(K,t 1 


A 

30 

00204 

334 


DO 40 J«3,NN 


A 

31 

00207 

344 


JMlaJ<*t 


A 

32 

00210 

354 


40 6(J,K)aG(2,K)4G(JHt,Kl 


A 

33 

00212 

364 


IF (NO) 50,70,50 


A 

34 

00215 

374 


60 00 60 Ja 1 ,N0 


A 

35 

00220 

384 


60 0( J4NN,K)a«G( J4l ,E)4T(K1 


A 

36 

00222 

394 


70 CONTINUE 


A 

37 

00222 

404 

c 

part Ml COMPUTE DOT PRODUCTS 


A 

38 

00223 

414 


DO 90 Jal ,N 


A 

39 

00226 

424 

' 

DO 80 lal , J 


A 

40 

00231 

434 


GDG( 1 ,J)«G( I ,K)4G( J,K)4G06( I ,U) 


A 

41 

00232 

444 


80 G0G(J,naG0G(t,Jl 


A 

42 

00234 

454 


90 YOG( J)aY(K)4GU,K)4YDG(Jl 


A 

43 

00236 

464 


Y0Yay(K|4Y(K)4Y0Y 


A 

44 

00237 

474 


lOO YOfJaYOCl.TT 


A 

45 

00237 

484 

c 

PART 2 MAXIMUM H(l) 


A 

46 

00241 

494 


XBL»1 


A 

47 

00242 

504 


HM«0« 


A 

48 

00243 

514 


DO 120 laK8L,N 


A 

49 

00246 

524 


B( 1 , 1 )4GD6( 1 ,I ) 


A 

50 

00247 

534 


A ( 1 , 1 )aY06( l)/eu , 1 1 


A 

51 

00250 

544 


D00«A( 1 , 1 )4Y0G( 1 1 


A 

52 

00251 

554 


H«0D0 


A 

53 

00252 

564 


IF (ABS(H)«ABS(HM11 120,120,110 


A 

54 

00255 

574 


1 10 BMF( 1 )aB( 1 • 1 ) 


A 

55 

00256 

564 


CA( 1 )aA( 1 , 1 ) 


A 

56 

00257 

594 


1M(1 )■! 


A 

57 
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00260 

60* 


DDDMaODO 

0026 i 

6 1 • 


HMaH 

00262 

62« 

120 

CONTINUE 

00268 

63* 


6(1,1 )«RhF( 1 1 

00265 

68* 


A ( 1 • 1 )«CA( ! ) 

00266 

65* 


000»1 •-OODM/YOY 

00267 

66* 


ODD^ABS ( DOO 1 

00270 

67« 


RMSaSQRT(000/P0INT» 

0027 1 

68* 


PRINT 820, K8L , IM(! ) ,CA( 1 1 ,0D0,RMS,HM 

00271 

69* 

C 

PART 3 maximum H(21 

00301 

70» 


KBLa2 

00302 

71* 


HMaO, 

00303 

72* 


00 150 Ial,N 

00306 

73* 


IF U-IM(l) 1 130,150,130 

0031 1 

78* 

130 

CONTINUE 

00312 

75* 


MaIM( 1 1 

00313 

76* 


8(2,1 )"6D6( I ,H) 

00318 

77* 


B( I ,2)aB(2, 1 )/B( 1 , 1 ) 

00315 

78* 


B ( 2i2 ) "GOG ( 1 , 1 ) -B ( 2 , M 8B ( 1 ,2) 

00316 

79« 


A(2,2)aYDG(l)-B(2,l)8A(t,M 

00317 

80* 


A(2,2)«A(2, 21/8(2,2) 

0O32O 

81* 


A (2, 1 )>A ( 1 , M-A (2,2)8B( 1 ,2) 

00321 

62* 


ODD- A (2, 1 )8Y06(M)8A (2,2) 8YDG( 1 ) 

00322 

83« 


HaOOD 

00323 

88* 


IF (ABS(H)-ABS(HM) ) 150,150,180 

00326 

85* 

I 80 

BM8( I )»B(2, 1 ) 

00327 

66* 


BMF( 1 )«B( 1 ,2) 

00330 

87* 


BMF(2)aB(2,2) 

00331 

88* 


IM(2)al 

00332 

69* 


CA(2)aA(2,2) 

00333 

90« 


CAM )«A(2, 1 ) 

00338 

91* 


HMaH 

00335 

92* 


OODMaDOO 

00336 

93# 

150 

CONTINUE 

00380 

98* 


B( 2, 1 }"8mB( 1 ) 

00381 

95* 


B( I ,2)>BMF( 1 ) 

00382 

96* 


B(2i2)«BMF(2) 

00383 

97* 


A(2,2)»CA(2) 

00388 

98« 


A (2i'l }"CA M ) 

00385 

99# 


DOOal ,a0D0M/Y0Y 

00386 

100* 


DDOeABS(OOO) 

00387 

101* 


RMS*S0RT(0D0/P0INT) 

00350 

102* 


PRINT 830, KBL,!H(2) 

00358 

103> 


DO 160 1*1 ,2 

00357 

1088 

160 

PRINT 850, I ,CA ( M 

00368 

105* 


PRINT 860, DD0,RMS,HM 

00371 

1068 


IF (NP-KBL) 180,180»!70 

00378 

1078 

170 

PRINT 870 

00378 

1 088 

C 

PART 88898 8A 88888 

00376 

1098 

180 

L*tM( 1 ) 

00377 

1 108 


DO 370 KBL«3,NT 

00802 

Ml* 


HMaO, 

00803 

1 128 


KBLMlaKBL-I 

00808 

1138 


DO 270 I«1 ,N 

00807 

1 188 


DO 190 Jal ,K8LM| 

00812 

1158 


IF (I-|M(JM 190,270,190 


A 58 
A 59 
A 60 
A 61 
A 62 
A 63 
A 6*1 
A 65 
A 66 
A 67 
A 66 
A 69 
A 70 
A 71 
A 72 
A 73 
A 7** 
A 75 
A 76 
A 77 
A 78 
A 79 
A 60 
A 61 
A 82 
A 83 
A as 
A 85 
A 66 
A 87 
A 86 
A 69 
A 90 
A 91 
A 92 
A 93 
A 9‘f 
A 95 
A 96 
A 97 
A 98 
A 99 
A 100 
A 101 
A 102 
A 103 
A 108 
A 105 
A 106 
A 107 
A 106 
A 109 
A 110 
A 111 
A 112 
A 113 
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001IS ll*» 1*0 CONTINUE 

00N17 IIT* B (KBL I I I *GD6 ( I |L I 

00*20 I|8* B ( t I KBL ) «0 I KBL I I I /B ( 1 I I J 

00*21 1|9« DO 210 J«2,KBLMI 

00*2* 120* 

00*25 121* 

00*25 122* fl ( KBL f U) *GOG I I (N I 

00*27 123* DO 200 K*liJNl 

Q0N32 !2** 200 B ( KBL i J I "B ( KBL i J I I K i J ) *B ( KBL i K I 

00*3* 125* 210 B ( J |KBL ) «B I KBL I J I /B ( J I J I 

00*35 I25« B ( KBL |KBL ) »GDG ( I 1 1 ) 

00*37 127» »<KBL>KBL)*T0G( I I 

00**0 128» DO 220 K«1,KBLN1 

00**3 I2*« BIKBL|K8L)-B(KBL,KBL»-BIK,KBL)»B(KBL,K) 

00*** 130» 220 AtKBLiKBL)»5(KBLiKBL)-»IK,K)*B(KBL|K) 

0O**5 131* A(KBL>KBL)>5IKBL.KBL)/BIKBLiKBLI 

00**7 I32» DO 230 J«liKBLMl 

00*52 133* KBMJ«KBL»J 

00*53 13*« MKBL,KBMJ)«A(KBNJ,KBMJI 

00*5* 135* DO 230 K»1|J 

00*57 135* KMJPK»K8MJ*K 

00*50 137* 230 A ( KBL , KBM J I • A ( KBL , KBM J I - A ( KBL , K« JBK ) *8 I KBM J , KM JPK 1 

00*53 138* OOP«Ot 

00*5* 139* DO 2*0 J*I|KBLM1 


00H67 

I10* 


MatH( J1 

00<470 

1«*1* 

2*40 

ODDbDOD^A (KBL f J)*V06(M) 

00*»72 

1<*2« 


D00"00D«A(KeL|KBU)»YD6M) 

00*»73 

1**3* 


H»DDO 

00H7M 

1 ‘I*** 


IF (ABS(H)-ABS(HHn 270»270|250 

60‘»77 

1 ^5• 

250 

IM(KBL)"I 

oo&oo 

1 1|5« 


HHaH 

OOSOl 

l**7* 


DOOMsDDD 

00S02 

i«ie* 


DO 250 Jal tKBLMl 

0060S 



BHRU)aB(KBL,J) 

OOS06 

I50» 


BHFl JlaB( JfICBL) 

00507 

151* 

250 

CA(J)«A(KBL|J) 

0051 1 

152» 


BMFlKBUaB(KBL,KBU 

005t2 

153* 


CA1KBU"A(KBL|KBL) 

00513 

15*»* 

270 

continue 

00515 

155« 


00 260 Jm } (KBLHI 

00520 

155* 


BfKBLiJl"BM8( J) 

00521 

isr* 


BU.KBU*BHF( J) 

00522 

15B» 

2B0 

ATFBti JlaCAf J) 

0052** 

!S9« 


B(KBL|KBLl*BNF(KBt) 

00525 

150* 


A1KBL|KBL)«CA1KBL) 

00525 

151* 


000«1 ••DDOM/YOY 

00527 

152* 


PRINT *f**0i KBLiIHlKBLI 

00533 

152* 


DO 290 1 1«1 |KBL 

00535 

155» 

290 

PRINT *ISOi I I |CA (I p 

b05**3 

155* 


ODDbABS(DOD) 

005«l*t 

155* 


RHS«SQRT(OOD/POlNT 1 

005«I5 

I57» 


PRINT **50» DOO,RHS«HH 

00552 

15B* 


IF {NP«KBL> 310«3t0i300 

00555 

159* 

300 

PRINT *»70 

00557 

170* 

310 

IF nc6L*NP) 370»320|320 

00552 

I7|» 

320 

DODaOt 


A II* 

A 115 
A 115 
A 117 
A 118 
A 119 
A 120 
A 121 
A 122 
A 123 
A 12* 

A 125 
A 125 
A 127 
A 128 
A 129 
A 130 
A 131 
A 132 
A 133 
A 13* 

A 135 
A 135 
A 137 
A 138 
A 139 
A 1*0 
A 1*1 
A 1*2 
A 1*3 
A I** 

A 1*5 
A 1*5 
A 1*7 
A 1*8 
A 1*9 
A ISO 
A 151 
A 152 
A 153 
A 15* 
A 155 
A 155 
A 157 
A 15B 
A 159 
A 150 
A 151 
A 152 
A 153 
A 15* 
A 155 
A 155 
A 157 
A 158 
A 159 
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00SA3 

172« 

00 360 JaltlPOINT 

A 

l»0 

00566 

173* 

ENCODE I560,LABELXM 01 )KBL 

A 

in 

00571 

176* 

ANUMaQ, 

A 

17* 

00572 

175» 

DENOMat • 

A 

I7J 

00573 

176* 

DO 350 (at |KBL 

A 

17^1 

00576 

177* 

NafHf 1 1 

A 

175 

00577 

178* 

ir fN«NNl 330,330»360 

A 

17* 

00602 

177* 

330 ANUHaANUMaCAf na6(K»J) 

A 

177 

00603 

180* 

60 TO 350 

A 

178 

0060*1 

181* 

360 DENOMaDCNONaCAl I 1 aG { N.NNa 1 » J ) 

A 

|7» 

00605 

182* 

350 CONTINUE 

A 

180 

00607 

183* 

CYvANUN/OENOH 

A 

181 

00610 

I86« 

CY«CT«TT 

A 

182 

0061 1 

185* 

OlfEaYC Jl-CY 

A 

183 

00612 

186* 

6oir«oirr/YU) 

A 

|8<l 

00613 

187* 

OOOaOODaDIEEaOirF 

A 

185 

0O616 

188a 

PRINT 680i T(J)|X(J,|}»CV,Oirr«PDir«DENOH 

A 

18* 

0062a 

189a 

LPLOTYUlaVIJI 

A 

187 

0O625 

|90a 

LPLOTX(J|al(J»l 1 

A 

188 

00626 

191* 

LPtOTZIJI-CY 

A 

18* 

00627 

192a 

360 CONTINUE 

A 

1*0 

00631 

193a 

CALL 8UIK3V < • 1 | 35 *L ABELX »LABCL Y • • I PO I NT » LPLOTX t LPLOT Y ) 

A 

1*1 

00632 

1 96a 

call 0UIK3V (0,6Q,LABELX,LABELY««IP0INTiLPL0TX|LPL0TZ) 

A 

1*2 

00633 

195a 

ODOaABS(ODO) 

'A 

1*3 

00636 

196a 

RMSaSQRT (DOO/POINT I 

A 

1*1 

00635 

197a 

ERRlaOOO/YOY 

A 

1*5 

00636 

I98a 

ERRORaSORTf ERR 1 /POINT) 

A 

1** 

00637 

199a 

PRINT 690i 00O|RNS«Y0Y*CRRt (ERROR 

A 

1*7 

00666 

200a 

370 CONTINUE 

A 

1*8 

00650 

20la 

NCaNCal 

A 

1** 

00651 

202a 

00 TO 10 

A 

200 

00652 

203a 

375 call ENOJOB 

A 

201 

00653 

206a 

STOP 

A 

202 

00653 

206a 

C 

A 

203 

00656 

206a 

380 FORHAT IIHI ) 

A 

20* 

00655 

207a 

390 FORNAT n6(0l6«8(I6«l6(I6(0l6«8( 12) 

A 

205 

00656 

208a 

600 FORNAT < 12A6) 

A 

20* 

00657 

209a 

610 format (2015.8) 

A 

207 

00660 

210a 

620 FORNAT < 1 HO « 30HAPPROX 1 NAT I NO FUNCTION NUMBER I2|6H IS G()2»1H),/// 

A 

208 

00660 

21 la 

1 1 tH0|7HCA(l ) •e23«|6(/(lH0(7HD0D «E23 » 1 6 » 6X » 6HRNS a£23 • 1 6 » 6X i SH 

A 

20* 

00660 

2l2a 

2HM aC23«|6t///t65(2Haa)»/////) 

A 

210 

00661 

2l3a 

630 format < lH0(30HAPPROXtNATtNG FUNCTION NUMBER 12»6H IS 6I12»1H)»/// 

A 

211 

00661 

216a 

!1 

A 

212 

00662 

21Sa 

660 format UH0«30HAPPR0X|MATING FUNCTION NUMBER I2»6H IS 6(12(1H),/// 

A 

213 

00662 

216a 

1 1 

A 

2I* 

00663 

217a 

650 format UH0*3MCA( (I2*3H) «(E23tU) 

A 

215 

00666 

218a 

660 format UHOiBHOOD aE23* 1 6 • 6X »7HRMS ' *E23 » 1 6 • 6X (5HHM »E23»|6(/ 

A 

21 * 

00666 

2|9a 

I///I 

A 

217 

00665 

220a 

670 format USfZHaa) •/////) 

A 

218 

00666 

221a 

680 format (IHOfZX.SHY aC 10*6 ,2X (5HX 1 •£ 1 0 • 6 » 2X »5HC V •£ 1 6 » 8 |2X |5HD 

A 

21 * 

00666 

222a 

1 irFaC16«8«2X(SHP0IFaE|6«8(2X(7HDEN0M •C)6»8) 

A 

220 

00667 

223* 

690 format nHOi///|2X(5HODO •£ 1 5* 8 • 2X , 5HRMS «E 1 S • 8 *2X « 5HYDY •E15«8«2X 

A 

221 

00667 

226a 

l•6HeRRl •eiS«6|2X»7MCRROR aC 1 5 *8 » /// »65 1 2Haa ) , /// ) 

A 

222 

00670 

22Sa 

500 format (2S(2Haa| »7X(|3HCASE NUMBER 1 1 2 * 1 H ) , 7X * 25 1 2Haa ) , //// ) 

A 

223 

00671 

226a 

510 format MH0f66HTHE NUMBER OF DATA POINTS USED |N THIS FIT IS |6i// 

A 

22 * 

00671 

227a 

1iII«66HTHC maximum NUMBER OF APPROXIMATING FUNCTIONS IN THE NUMERA 

A 

225 
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